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ABSTRACT 

We develop a general method to fit the planetary distribution function (PLDF) to exoplanet survey 
data. This maximum likelihood method accommodates more than one planet per star and any number 
of planet or target star properties. Application to Kepler data relies on estimates of the efficiency 
of discovering transits around Solar type stars by Howard et al. (2011). These estimates are shown 
to agree with theoretical predictions for an ideal transit survey. Using announced Kepler planet 
candidates, we fit the PLDF as a joint powerlaw in planet radius, down to 0.5-R®, and orbital period, 
up to 50 days. The estimated number of planets per star in this sample is ~ 0.7 — 1.4, where the 
broad range covers systematic uncertainties in the detection efficiency. To analyze trends in the PLDF 
we consider four planet samples, divided between shorter and longer periods at 7 days and between 
large and small radii at 3 i?®. At longer periods, the size distribution of the small planets, with index 
a ~ —1.2 ± 0.2 steepens to a ~ —2.0 ± 0.2 for the larger planet sample. For shorter periods, the 
opposite is seen: smaller planets follow a steep powerlaw, a ~ — 1.9 ± 0.2 that is much shallower, 
a ~ —0.7 ± 0.2 at large radii. The observed deficit of intermediate-sized planets at the shortest 
periods may arise from the evaporation and sublimation of Neptune and Saturn-like planets. If the 
trend and explanation hold, it would be spectacular observational confirmation of the core accretion 
and migration hypotheses, and allow refinement of these theories. 

Subject headings: Methods: statistical — Planetary Systems — Planets and satellites: detection — 
Planets and satellites: dynamical evolution and stability — Planets and satellites: 
formation — Stars: statistics 



1. INTRODUCTION 

Individual exoplanet discoveries highlight the extraor- 
dinary diversity of worlds in the Solar neighborhood 
(Mayor & Queloz 1995; Marcy & Butler 1996; Charbon- 
neau et al. 2009; Lissauer et al. 2011a). For an accurate 
census of the planet population, the statistical analy- 
sis of large samples of exoplanets is required (Cumming 
et al. 2008; Howard et al. 2010). Trends revealed by the 
data provide powerful constraints on dynamical theories 
of planet formation (Goldreich et al. 2004; Kenyon & 
Bromley 2006; Youdin 2010). 

The correlation of giant planets with host star metal- 
licity is perhaps the most interesting trend revealed by 
radial velocity surveys (Gonzalez 1997; Fischer & Valenti 
2005; Johnson et al. 2010). This trend has been a pow- 
erful guide to identifying the mechanisms responsible for 
planetesimal formation (Youdin & Shu 2002; Youdin & 
Goodman 2005; Johansen et al. 2009), see Chiang & 
Youdin (2010) for a review. 

The Kepler transit survey is currently revolutionizing 
the field of exolanets from space (Borucki et al. 2010). 
Borucki et al. (2011, hereafter BKepll) released 1,235 
partially vetted planet "candidates." Morton & Johnson 
(2011) estimate the false positive rate as < 5% around 
the brighter stars that are considered here. For brevity 
and statistical purposes, we mostly refer to the candi- 
dates as "planets," though significant followup work re- 
mains with only 18 candidates currently confirmed. 

Howard et al. (2011, hereafter HKepll) presented a 
detailed statistical analysis of the ~ 500 BKepll planets 
around ~ 60, 000 bright, Solar-type stars. Accounting for 
detection efficiencies, HKepll report planet occurrence 



rates of ~ 0.2 planets per star for radii > 2i?® and orbital 
periods < 50 days. The purpose of this paper is to apply 
different statistical methods to the Kepler dataset, mak- 
ing use of the estimates of detection efficiencies — whose 
importance is on par with the detections themselves — 
provided by HKepll. 

Our method is based on the Tabachnik & Tremaine 
(2002, hereafter TT02) technique to analyze the mass 
and period distributions of planets from radial velocity 
surveys. TT02 emphasized the importance of determin- 
ing planet mass and period distributions simultaneously 
A simultaneous fit is crucial because radial velocity de- 
tection thresholds depend on both the mass (times the 
sine of inclination) and period. Transit surveys have a 
similar, though weaker, coupling between planet radius 
and orbital period in the efficiency of discovering transits. 

A primary advantage of the TT02 technique is that it 
naturally accommodates more that one planet per star. 
Planet hosting is not treated as a binary proposition. 
This distinction is not just a technical nicety, because the 
number of planets per star is order unity at least. We 
show that Kepler data already imply about one planet 
per Solar type star within ~ 0.25 AU and with a radius 
bigger than 0.5i?®. This number is certain to grow with 
time as longer periods and smaller sizes are probed. 

High multiplicity rates in the Kepler sample (Latham 
et al. 2011) means that the fraction of stars that host 
a planetary system (at most one, of course) can be 
significantly smaller than the number of planets per 
star (NPPS). We treat the planet population (includ- 
ing known multiples) as a whole, and ignore multiplicity 
issues. Ragozzine & Holman (2010) discuss the value of 
multi-transiting systems in constraining mutual inclina- 
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tions. Using these constraints, Kepler data show that 
the number of planets per planetary system is at least 2 
— 3 (Lissauer et al. 2011b). 

This paper is organized as follows. Section 2 describes 
the detection efficiencies relevant to the Kepler transit 
survey, including an analytic fit to discovery efficiencies 
reported by HKepll in §2.2. Section 3 describes our 
method for analyzing cxoplanct survey data. Section 
4 applies the method by fitting Kepler data to a joint 
powerlaw in radius and period. Differences between the 
shorter and longer period planets (in §4.2) and also be- 
tween the smaller and larger planets (in §4.3) are ana- 
lyzed. Section 5 gives non-parametric estimates of planet 
occurrence, and comments speculatively on the conse- 
quences of rising period distributions. We conclude with 
a summary and discussion in §6. Appendix A gives an 
analytic solution for powerlaw fits to data from an ideal- 
ized transit survey. 

2. SELECTION EFFECTS FOR TRANSIT SURVEYS 

A robust statistical analysis must account for relevant 
selection effects. We quantify selection effects as detec- 
tion efficiencies, r), which give the ratio of detections to 
actual planets. Individual efficiencies multiply to give 
the net detection efficiency. 

The three main selection effects for transit surveys are: 
(i) the transit probability for a planetary orbit to cross 
our line of sight to the star, n tI < 1; (ii) the discov- 
ery efficiency of the survey in finding transiting planets, 
?7disc < 1; and (iii) the rate of false positive events that 
mimic a planet transit, n p < 1. For false positives, the 
relevant r/f p = 1/(1 — r rp ) > 1, the only efficiency that 
can exceed unity. As mentioned already, the false posi- 
tive rate is estimated to be low enough that we ignore it 
for simplicity. 

For most analyses, an average detection efficiency is in- 
sufficient. The dependence of the efficiencies on relevant 
properties of planets and target stars is required. We 
consider how the efficiencies vary with planet radius, R, 
and orbital period, P. See §3.3 for a general discussion 
of which parameters should be included. 

The left panel of Fig. 1 shows the Kepler discovery 
efficiency, calculated as described below (in section 2.2). 
The discovery efficiency is nearly unity over a signifiant 
range of planetary radii and periods, thanks to the high 
photometric precision of Kepler. There is a sharp drop 
in discovery efficiency for sufficiently small planets whose 
transit depths compete with noise in the photometric 
data. At shorter periods, the signal from smaller planets 
can rise above the noise because more transits are seen. 
The right panel of Fig. 1 shows how the geometic transit 
probability reduces the probability of finding long period 
planets. 

To connect with the extensive work on radial velocity 
(RV) surveys, we note that RV upper limits can be ex- 
pressed as a detection efficiency for a given star. Most 
simply, r/ = 1 for planet parameters above the detection 
threshold and rj = otherwise. This step function can 
be smoothed to account for confidence intervals. 

2.1. Transit Probability 

The probability that a planet, on a randomly oriented 
circular orbit of semi-major axis a, transits it's host star, 



with radius R* and mean density p*, is 

for R <C i?*. Eccentricity tends to increase the de- 
tectability of planets at fixed a (Burke 2008). 1 For 
the mean eccentricities ~ 0.1 — 0.25 implied by Kepler 
transit durations (Moor head et al. 2011), corrections are 
< 10% and ignored here. 
For simplicity, we approximate the mean stellar density 

(or more specifically the mean 1 ^ 3 ) as Solar. The den- 
sities of planet hosts can be estimated from precise light 
curve parameters, though corrections for the eccentricity 
apply (Tingley et al. 2011). Furthermore, planet hosts 
have a lower average density than the target star popu- 
lation, on general. This bias arises because lower density 
stars have higher transit probabilities. For our purposes, 
uncertainties in the average stellar density modestly af- 
fects the magnitude (but not the shape) of the inferred 
planet distribution. 

One might (incorrectly) expect multi-planet systems 
to require special values of r/tr- When mutual orbital 
inclinations are small, finding one planet does increase 
the odds that others will be found. However since the 
planetary systems as a whole are randomly oriented, low 
mutual inclinations also make it easier to miss an en- 
tire planetary system. The mutual inclination of planets 
within systems has no effect on the overall detection rate 
(aside from sampling noise, as always). 

2.2. The Kepler Discovery Efficiency 

Precisely quantifying the discovery efficiency of a sur- 
vey is quite difficult and is ongoing work for the Ke- 
pler mission. For our study we rely on the estimates of 
discovery efficiency in HKepll. HKepll quantified the 
discovery efficiency for a subsample of bright solar-type 
stars with relatively high r^isc- The sample consists of 
N* = 58, 041 stars that satisfy the effective temperature, 
surface gravity and Kepler magnitude cuts: 

4100 K < r cff < 6100 K (2a) 

4.0 < log 5 /(cm s~ 2 ) < 4.9 (2b) 

K P < 15 mag . (2c) 

Section 2.3 discusses the planet candidates in this sam- 
ple. 

The data for r;<jisc is presented in Figure 4 of HKepll. 
Here the discovery efficiency is reported for planets with 
0.68 < P/(days) < 50 and 1 < R/R m < 16 on a log- 
uniform grid. HKepll report rydisc-Y*, the number of 
stars around which a planet with that radius and period 
could be detected, in the bottom left of each cell. The 
related values of ?7disc are plotted in Fig. 2. 

For bins with no detected planets, HKepll do not 
report a discovery efficiency. This omission is not be- 
cause planet detections are required to estimate efficien- 
cies. Rather HKepll applied efficiencies to the detec- 
tions only, which differs from our approach of applying 

1 For an instructive explanation see http://oklo.org/2011/03/ 
23/the-eccentricity-distribution/ 
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Fig. 1. — (Left): The Kepler mission's efficiency of discovering transiting planets around bright Solar-type stars vs. planetary radius 
and orbital period. The discovery efficiency is unity for large planets, but drops sharply below a planet size that increases gradually with 
orbital period. See text for details. (Right): The net detection efficiency combines the discovery efficiency (at left) with the geometric 
transit probability, which exerts a bias against detecting long period planets of all sizes. 



I 



10 4 
10 3 
10 2 
10 

1 

10" 1 
10" 2 



•Hh 



+ 
+ 
+ 

+ 
+ 



'•+. 



"■■+•+, 



16</? p < 16*2 1/2 " 
+;!|B*2 1/2 </? p < 16 
*"-8<# p <8*2 1/2 
+1 4*2 1/2 <7? n < 



'+. 



+. 



+■•.+ 4</? p <4*2 



8 

1/2 



•+. 



2*2 1/2 < R p < 4 



+ 2 < R D < 2*2 1/2 . 



+. 



2 1/2 < R p < 2 



1 <R P <2 



1/2" 



10 100 
Period [days] 



1000 





10 4 - 




10 3 - 


o 
(/) 


10 2 - 


I 






10 - 








1 - 




10" 1 - 




10" 2 




1 10 
Planet Radius [Earth Radii] 



Fig. 2. — The Kepler discovery efficiency plotted as a ratio, ??diac/(l — ?7disc) of detectable-to-non-detectable planets. Symbols correspond 
to the values reported in Figure 4 of HKepll, while the curves give the broken powerlaw fit of Equations (4) and (6). At left, the ratio 
is plotted vs. orbital period for a range of (binned) planetary radii. The same data are plotted vs. planet radius at right. The smooth 
behavior of the Kepler discovery efficiencies is evident. 



the efficiencies across all parameter space. The missing 
efficiencies in Figure 4 of HKepll can be readily obtained 
by interpolation. Most of the empty bins correspond to 
small P and large R where 77di sc ~ 1- Thus the lack of 
detections in these bins carries heavy statistical signifi- 
cance. 

Instead of merely interpolating the reported Kepler ef- 
ficiencies, we are motivated by their smooth variation to 
find an analytic fit. We find an excellent joint powerlaw 



fit, not to 77disc itself, but to the related 



Vdis 



1 



the ratio of discoverable to non-discoverable planets. 
Since r^sc ranges from zero to arbitrarily large values, 
it is simpler to fit with a powerlaw. The restriction to 
%isc < 1 is automatically satisfied. 
Fig. 2 shows that the broken powerlaw is an excellent 
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fit to reported discovery efficiencies. The powerlaw fit is 
/ R \ 6 ' 34 / P V 107 

below a break at 

/ P \ 017 

R < Rb = 2.2 . (5) 

The efficiency at the break is quite high, rydisc = 0.98 
(and remarkably constant, probably reflecting the high 
S/N threshold). Note that 7/disc drops to 50% (rdi sc = 1 
in Fig. 2) at R = 0.54i?b, rising roughly from 1 to 2 R® 
over the periods considered. When discovery efficiencies 
are small, 7/disc ~ ^disc *C 1, the powerlaw in Equation (4) 
applies directly to ?7disc- 

While the scalings for ?ydi S c > 0.98 are of little practical 
concern for our study, we report the fit above the break 
for completeness: 

/ R \ 2 - 74 / P \ -° Aeo 
^>^ b ) = 5.46(-j (-) . (6) 

This fit describes a small fraction of noisy stars, but does 
not affect our results because ?7disc ~ 1 in this regime. 

The relevant fit in Equation (4) agrees well with theo- 
retical estimates. For S/N limited detections, the com- 
bined efficiency should scale as 

?7disc77tr CX p- 5 / 3 i? 6 , (7) 

see e.g. Equation (7) of Gaudi (2007). The small R limit 
of Equation (4) gives r^isc^tr oc P~ 1 - 74 R 6 M , rather good 
agreement. 

The efficiencies summarized here will likely be im- 
proved by more detailed studies that include extraction 
of simulated transits from the actual Kepler data analy- 
sis pipeline. Accurate estimates of the Kepler discovery 
efficiency for all relevant parameters — especially stellar 
T e ff — are the most crucial ingredient for determining 
the frequency of Earth-like planets. 

2.3. The Planetary Sample 

Tables 1 and 2 of BKepll list the properties of 1,235 
planet candidates and their host stars. We must restrict 
our attention to the planet candidates around host stars 
with known discovery efficiencies, rydisc- Thus we only 
consider host stars within the Solar sample of HKepll, 
set by Equation (2). We further restrict attention to de- 
tections with signal-to-noise (S/N) > 10 and P < 50 
days, as these are also conditions for the validity of ?7disc- 
Applying these filters reduces the number of planet can- 
didates from 1,235 to 566. 

Our "full" planet sample considers 

0.5 days < P < 50 days (8a) 
0.5 i? ffi < R < 20 R® . (8b) 

Adding these radius and minimum period cuts only re- 
moves four planets (from 566 to 562). Instead of us- 
ing the BKepll reported values for R, we compute R 
from the reported values for the host star radius, i?*, 
and R/R* 2 

2 This step overcomes the rounding of the reported R to 0.1R® 



We now discuss relevant differences between our planet 
sample and that of HKepll. HKepll conservatively re- 
stricted attention to R > 2R®, because the discovery 
efficiency is high for these planets (except at longer peri- 
ods). We extend our analysis down to R = 0.5R®, which 
includes essentially all small planets in the Solar sample. 
This extension involves (a) trusting the HKepll reported 
efficiencies down to the 1-R© level to which they were re- 
ported, and (b) extrapolating the efficiencies down to 
0.5.R®. As shown above, both the smooth variation of 
the reported efficiencies and their agreement with the- 
ory give us confidence in making the extrapolation. We 
also report results for samples of larger planets that are 
unaffected by this extrapolation. 

HKepll defined S/N slightly differently than the pub- 
licly reported BKepll S/N values on which we rely. 
HKepll define S/N based on only one quarter of data 
(Q2, Kepler's first full quarter). By contrast BKepll 
report S/N values up to Q5, so that roughly four times 
more data has been collected. With perfect noise scaling, 
the BKepll S/N values should be roughly twice as high 
as those used by HKepll. Table 2 of HKepll presents 
noise scalings for 4 planets with ~ 2.5R® radii and ~ 40 
day periods. The ideal scaling roughly holds for three 
objects, but S/N saturated and only rose modestly for 
the fourth. It's unclear what the general noise scaling is, 
and how it varies with size and period. 

We proceed conservatively by performing all our anal- 
yses on both a S/N > 10 and a S/N > 20 planet sample 
(as defined by BKep 11). These cases cover the com- 
plete range of possibilities between saturated noise and 
perfect scaling. Fortunately the adopted S/N threshold 
has a modest effect on quantitative results (most impor- 
tantly, the number of planets per star) but no effect on 
the interesting trends we discuss. 

We briefly note a discrepancy between the planet 
counts reported by BKepll and HKepll. HKepll re- 
port a sample of 438 planet candidates (of the 1,235 re- 
leased by BKepll) that (a) have hosts that satisfy the 
stellar parameter cuts of Equation (2), (b) have planet 
parameters P < 50 days and 2i? ffi < R < 32R® and (c) 
satisfy the HKepll definition of S/N > 10. If we apply 
cuts (a) and (b) only to the BKepll tables, but allow all 
S/N, there are 378 candidates. That number increases 
to 393 using the reported R instead of calculating it as 
(i?/i?*) x Thus there are at least 60 (or 45) planets 
that HKepll say should appear in the BKepll tables, 
but they do not. 

Applying any S/N threshold to the BKepll tables can 
only make the discrepancy larger. We made no choices 
that would cause us to miss planets in the BKepll ta- 
bles. Cuts were applied inclusively (including any val- 
ues that fall on a boundary) and all the planets in any 
system were counted. We conclude that a discrepancy 
exists between the planet and/or stellar parameters used 
by HKepll and reported by BKepll. This discrepancy 
does not affect our results as long as the estimates of the 
discovery efficiencies in HKepll are sufficiently accurate. 

3. A GENERAL METHOD FOR EXOPLANET STATISTICS 



intervals, and is merely for convenience in defining boundaries in 
an assumption-free way. 
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This section describes a general method for estimat- 
ing the planetary distribution function (PLDF) from the 
results of a survey with quantified detection efficiencies. 
The purpose of developing this formalism is twofold. 

First, this method can investigate a PLDF that de- 
pends on properties of both the planet detections and 
the target stars. Further, the PLDF can take an arbi- 
trary functional form. While the exoplanet community 
is unlikely to agree on a single statistical methodology 
anytime soon, at least such an approach is possible. 

Second, this formalism allows the inclusion of all de- 
tected planets — including those in multi-planet systems. 
This inclusion is possible because we treat the PLDF as 
the average number of planets per star (NPPS) not the 
fraction of stars with planets (FSWP). Following TT02, 
our technique treats planet occurrence as a Poisson pro- 
cess, i.e. a series of independent random events. This 
technique naturally gives the NPPS. 

By contrast, many statistical analyses treat planet 
hosting as a binary proposition, i.e. a star either does 
or does not have a detected planet. The use of binomial 
statistics naturally gives the FSWP. Including multiple 
planets per system in this type of analysis is formally 
incorrect, though discrepancies are small if the multi- 
plicity rate is low. For more discussion of these points 
see Cumming et al. (2008), who rigorously analyze the 
FSWP in radial velocity surveys by only including the 
most detectable planet around any star. 

For transit surveys, the FSWP is more difficult to de- 
termine and is not addressed here. The main issue is 
applying the transit efficiency correction. Should the 
missed planets that do not transit be assigned evenly 
among all stars, or into high order multi-planet systems? 
As mentioned in the introduction, multi-transiting sys- 
tems are a powerful constraint, as are transit timing vari- 
ations and comparison to RV surveys (Ragozzine & Hol- 
man 2010). 

Before developing the general method for fits to a 
parameterized PLDF in Section 3.2, we consider non- 
parametric estimates of the NPPS in Section 3.1. Para- 
metric fits also give an estimate of the planet occurrence, 
whose quality depends on the appropriateness of the as- 
sumed functional form. The main reason to consider pa- 
rameterized fits is to identify and assess trends in the 
data — and if feeling bold, to extrapolate. Section 3.3 
discusses the consequences of fitting some parameters 
and ignoring others. 

3.1. N on- Parametric Estimates of Planet Occurrence 

For a survey of stars, we divide planet detections 
into bins indexed by £, with Ni planet detections in a 
given bin. If the average detection efficiency per bin is 
T)e, then the best estimate of the NPPS in each bin is 



fa 



(9) 



the number of detected planets per star divided by the 
efficiency. 

The total planet occurrence, y\ fa, depends on bin size 
when rji is not constant. For arbitrary small bins, there 
is only one planet per bin, and the efficiency is evaluated 
for the parameters of each planet and its host. This 
small bin limit is the only truly assumption-free way to 



estimate the planet fraction, and suggests that binning 
of data is never required. 

A maximum likelihood analysis using Poisson statistics 
can reproduce the estimates of fa given by Equation (9), 
and develop intuition for the general method. The num- 
ber of expected detections in all bins is 



N * ^2 Vefa ■ 



The likelihood function for a Poisson process is 



L 



Wimfa) 



CXp (-A^oxp) 



(10) 



(11) 



which represents the product of the probabilities of each 
individual detection (the term in square brackets) times 
the probability of finding no additional planets in any 
bin (the exponential). 

Any constant multiplicative factors can be ignored in 
the likelihood function, because they do not affect the 
maximization of likelihood. (Thus factorials do not ap- 
pear in Equation (11).) Here constant means indepen- 
dent of the PLDF, i.e. the fa values. Remarkably, this 
freedom allows us to ignore the efficiencies that charac- 
terize the Ni detections, but not in 7V cxp . The likelihood 
thus simplifies to 



CXp (-TVexp) 



(12) 



The efficiencies now appear only in iV exp . 

The maximum likelihood values of fa are the roots of 
dL/dfa = 0, which reproduces the expected result of 
Equation (9). It is easier (analytically and numerically) 
to maximize ln(L). 

A fundamental feature of the method is that efficien- 
cies are used to predict the number of expected planets 
and are not directly applied to the detections. 3 While 
this feature disappears when bin sizes are chosen to be 
arbitrarily small (as discussed above) , it becomes impor- 
tant for parametric fits. 

3.2. Fitting the Planetary Distribution Function 

We now describe a general method to fit a parame- 
terized PLDF to unbinned data. We consider a PLDF 
that depends on planet properties, x (an N x component 
vector) and stellar properties z (with N z components). 
Usually the elements of these vectors are logarithms of 
measured quantities such as the orbital period or stellar 
mass. 

The probability that a star with properties z has a 
planet with properties x that lie in a volume 4 dx of phase 
space is 

* = ^«- («) 
Defined this way, the integrated / represents the NPPS. 

3 Non-detections however can enter 7V exp via upper limits that 
define the efficiency, see §2. 

4 For brevity the exponent in the phase space volume is dropped, 
i.e. dx = dx Nlc here and throughout. We keep this exponent in 
the denominator of derivatives to be explicit that these are not 
gradients, as in Equation (13). 
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We express the differential distribution, 

d J^^C 9 { X ,z; a) , (14) 

as an amplitude, C, times a shape function, g. The shape 
parameters, ot, can describe the behavior of individual 
planetary or stellar properties as well as any correlations. 
In the simplest case, the values of a are powerlaw expo- 
nents. 

The total number of planets around stars (indexed 
by j) is 



N, 

N to t = / g(x,z J ;a)dx 

j J 

7 



(15a) 



= N^C J F*(z)g(x, z:a)dxdz . (15b) 

where integration covers a specified volume of the phase 
space of x. In Equation (15b) the sum over stars is re- 
placed by an integral over the stellar distribution func- 
tion F*(z). We proceed with the integral notation to 
minimize indices, but direct summation over all stars is 
always possible (and preferred when feasible). 
Consider a survey with a net detection efficiency 



V(x,z) = Y[r] k (x,z) 



(16) 



that is a product of efficiencies from individual selection 
effects (indexed by k). The number of expected planet 
detections is 



N< 



cxp 



N*C J r}(x, z)T*(z)g(x, z; ct)dxdz , (17a) 

(17b) 



N*CF 



The shape integral, F, weights the shape function over 
both the stellar distribution and the efficiencies. 

3.2.1. Maximizing Likelihood 

Consider a survey that detects N p \ planets around N* 
stars. We define the likelihood of the data analogously 
to Equation (12) as 



exp(-7V oxp ) , 



(18) 



which represents the probabilities of each individual de- 
tection, labelled by i, times the probability that no other 
planets were detected, given by the exponential factor. 

Applying Equations (13) and (14), we again eliminate 
unnecessary constants, here the phase space volume, to 
define the likelihood function as 



C Npl Y[g{xi,z j{i y,a) 



cxp (-N cxp ) , (19) 



where refers to the star that hosts planet detection i. 
We let gi represent the shape function evaluated for the 
properties of detection i and express the log likelihood 

as 

In(L) = N pl ln(C) + £ ln( 5i ) - iV cxp . (20) 

i=l 



The best fit normalization, C, is found by maximizing 
the likelihood as the root of d In L/dC = 0: 



C* = 



N, 



N*F 



(21) 



using Equation (17b). 

Following TT02, we use this constraint to eliminate C 
from the likelihood: 



JVpi 

ln(L) = -7V pl ln(F)+^ln(. 9 

i 



(22) 



iV, 



pi 



In 



where the constant term in curly brackets can again be 
ignored. 

The best fit values of a. maximize L as the roots of 
dL/da = 0, which are 



dlnF 
da 



Afpi 



— E 

^Pi 2 r 



dln(gj) 
doc 



(23) 



These N a (the number of parameters in a) equations 
generally couple to each other and require numerical so- 
lution, but see §A. Equation (23) weights the shape func- 
tion by detection efficiencies on the left hand side and 
over individual detections on the right hand side. 

The basic steps in obtaining a best fit solution are 
as follows. First, select a form of the PLDF to fit the 
data and express in the form of Equation (14). Second, 
solve Equation (23) for the shape parameters as just de- 
scribed. Third, calculate the normalization via Equa- 
tion (21) with F evaluated using the best fit shape pa- 
rameters. The best fit solution is now complete and can 
be used (e.g.) to estimate the NPPS as N tot /N* using 
Equation (15). Fourth, estimate the errors on the fit via 
likelihood contours given by Equation (20), as described 
in §3.2.3 below. 

3.2.2. Interpretation 

This formalism has a remarkably straightforward in- 
terpretation. The normalization C matches the numbers 
of detected and expected planets, N p \ — N eKp , as seen by 
comparing Equations (17b) and (21). 

The shape parameters similarly match the expected 
and detected averages of planet and host star observ- 
ables. Consider the case of powerlaw distributions with 
a shape function 



exp (a p i • x + a* • z) 



(24) 



Here x and z are the logarithms of quantities being fit to 
a powerlaw, and a = {a p i,a*} distinguishes the pow- 
erlaw exponents for planetary and stellar parameters. 
With this powerlaw shape function Equation (23) gives 



(x) F = {x) obs , (z) F = {z) ohs 



(25) 



where 



{{x) F ,(z) F } = J {x,z}r]F*gdxdz (26a) 



JVp, 



{(a;) obs , (z)obs} = -jTj-^ixuZi} . 

Jv pi i 



(26b) 
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Non-powerlaw shape functions can investigate more de- 
tailed properties of the detections than simply the aver- 
age value of observables. 

3.2.3. Errors 

To quantify the uncertainty in the shape parameters, 
we compare contours of the likehood L to a Gaussian 
distribution. The maximum likehood, L max , is given by 
Equation (22) with the best fit a. The la uncertain- 
ties are defined by the ln(L) = ln(L max ) — 1/2 contour. 
The general na contours follow ln(L) = ln(L max ) — n 2 /2. 
The uncertainty in the normalization is just the range of 
values taken by C within the error ellipse of the a. 

The ID errors on an individual shape parameter (a 
component of a) are usually given assuming other pa- 
rameters are held at their best fit values. When there is 
strong covariance, and a very elongated error ellipse (not 
the case in this study), ID errors are not very meaning- 
ful. 

The errors on a. describe the precision of a fit, not the 
quality. A large planet sample precisely defines the aver- 
age powerlaw even if the actual PLDF deviates strongly 
from a powerlaw. The quality of fit can be determined 
by Kolmogov-Smirnoff (K-S) tests. 

This derivation ignores uncertainties in the detection 
efficiencies and the planetary and stellar parameters. 
Here, we describe how to include these errors but do 
not include them in our analysis. Uncertainties asso- 
ciated with an individual detection, i, affect the value 
of the shape function, gi, that appear in the likelihood 
analysis. To include these uncertainties, gi should be a 
weighted integral over the uncertainties 

9i = J <j)i{x' -x u z' - z j{ i))g{x',z';a)dx'dz' (27) 

where & is an appropriately normalized Gaussian distri- 
bution (for instance) centered on the best fit values of 
the planet detection and its host star. Ignoring uncer- 
tainties is equivalent to setting the error distributions, 
<pi, to 5- functions. 

Uncertainties in the target star properties affect the 
number of expected planets via the shape integral F. 
Errors can be included similarly to Equation (27) as 

1 n, 

F= J^Y,J -ZjHx,z')g(x,z')dxdz' (28a) 

* 3 

w J 4>{z' - z)T* (z)r)(x, z')g(x,z')dxdz'dz. (28b) 

The two forms above show how errors could be applied 
to each star j, or (more approximately) how average er- 
rors could be assigned to the stellar distribution function. 
The integration over both z and z' in Equation (28b) is 
potentially confusing. The integration over z is over the 
stellar distribution and replaces the sum over individual 
stars in Equation (28a), while the integral over z' covers 
the range of uncertainties allowed by the error distribu- 
tions cj). 

Estimates of the total number of planets should also 
include target star uncertainties, by similarly modifying 
Equation (15) (i.e. setting r] — >• 1 in Equation (28)). 

Systematic uncertainties in the detection efficiencies 
are more difficult to quantify. Assigning detection effi- 
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ciencies are the most important (and dangerous) part of 
any statistical analysis. It is safest to perform multiple 
analyses that cover a range of (hopefully reasonable and 
nearly complete) possibilities for detection efficiencies. 
That approach is the one taken here. 

3.3. Which Parameters to Study 

In any analysis one must choose to study some subset 
of planetary and stellar parameters and to ignore others. 
For a robust fit, the PLDF should include all parameters 
that cause the detection efficiency to change significantly 
across the parameter space studied. Ignoring a param- 
eter in the PLDF amounts to averaging the PLDF over 
the relevant parameter. However as expressed in Equa- 
tion (17), the correct way to predict the yield of a survey 
is to weight the PLDF by the detection efficiencies before 
averaging. 

A safer way to ignore a parameter (relevant to the de- 
tection efficiency) is to first fit it and then marginalize 
over it. With this caveat made explicit we now show how 
the general formalism treats the case of a PLDF that de- 
pends only on planetary parameters, x, or only on stellar 
parameters, z, (most famously stellar metallicity) . 

3.3.1. Planetary Parameters Only 

For a shape function, g(x;at), that is independent of 
stellar properties, the shape integral simplifies to 

F = J rj4x)g(x;a)dx. (29) 

where the stellar-averaged efficiencies are 

r]*(x) = J rj(x,z)J 7 ^(z)dz. (30) 

The fit procedure is simplified in this case. The stellar 
averaging is done a single time, and not for every choice 
of a. during optimization. 

This case is the one most relevant to this work. Specif- 
ically in Section 4 we consider a shape function 



g = cxp(ax + By) = 




This powerlaw distribution in planetary radius and or- 
bital period, identifies the planetary parameters as 

x = {x,y} = {ln(R/R )MP/Po)} (32) 

with the help of a reference radius R and period P , and 
defines the shape parameters a = {a, 6} as powerlaws. 
Appendix A gives an analytic solution for these best fit 
powerlaw exponents for an ideal transit survey with unity 
discovery efficiency and no (unvetted) false positives. 

3.3.2. Stellar Parameters Only 

We now consider fitting a PLDF that depends on stel- 
lar parameters only. As cautioned at the beginning of 
this section, it would be a bad idea to do this for a tran- 
sit survey since the detection efficiencies always vary with 
orbital period, and also vary with planetary radius un- 
less small planets with a discovery efficiency below unity 
are ignored. We consider this case for other survey types 
and for completeness. 
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Since the PLDF is independent of x by assumption, 
we can integrate the differential distribution of Equa- 
tion (14) over a volume, Vx, of planetary parameter 
space. Defining G = gVx gives 

f{z)=CG{z;cx*), (33) 

which represents the NPPS with characteristics z. The 
shape integral simplifies slightly to 

N cxp = N*C [(ri(z))F*{z)G(z,a*)dz (34) 

(r ] (z)}=fr](x,z)dx/V x , (35) 

by averaging the efficiencies over planetary parameters. 

The best fit solutions are still given by Equations (21) 
and (23), and \ngi can be replaced with lnG^ (no ap- 
proximation needed) to remove all references to the dif- 
ferential distribution. 

4. POWERLAW FITS TO KEPLER DATA 

We now use the Kepler planet candidates announced 
by BKcpll to constrain the underlying PLDF (planetary 
distribution function) as a powerlaw, 

d\nR3\nP =C (£) GO ' (36) 

in planet radius, R, and orbital period P. Equation (36) 
represents the NPPS (number of planets per star) per 
logarithmic interval in R and P. The actual PLDF 
can deviate from a powerlaw (and the data show that 
it does). Nevertheless, powerlaws are a common and 
useful approximation that captures basic trends in the 
data. Identifying deviations from powerlaw behavior re- 
veals the relevant scales of planet formation and migra- 
tion, and we show that Kepler data identify these trends 
extraordinarily well. 

We study the full planet sample around Solar type stars 
in §4.1. In §4.2, we compare the distributions of shorter 
and longer period planets. Further comparing results for 
larger and smaller planets in §4.3 shows that the sample 
has a deficit of intermediate-sized planets at the short- 
est orbital periods. For all these correlations, the S/N 
threshold has limited impact (§4.4). 

4.1. The Full Sample 

Our most complete sample of planets covers 0.5 < 
P/day < 50 and 0.5 < R/R e < 20. (See §2.3 for de- 
tails.) Fig. 3 plots planet counts of this full sample. The 
planet counts are binned by R (P) in the left (right, re- 
spectively) panel. Raw counts are divided by the number 
of stars surveyed and the logarithmic bin size. 5 Normal- 
ized this way, the data measure the PLDF as weighted 
by the detection efficiencies, i], 

where brackets indicate marginalization over P and R, 
respectively. Aside from sampling noise, these values are 
independent of the survey size or choice of bin width. 

5 E.g. divided by ln(ilj_|_i/.Rj) for a bin between Ri and 
and similarly for period bins. 



The best fit PLDF is "projected" into observational 
space in Fig. 3. This projection requires weighting by the 
detection efficiencies of transits, 7/t r , and discovery, fjdiso 
before marginalizing over P (or R), as in Equation (37). 
The comparison of the fits to histograms is somewhat 
misleading. The actual analysis uses unbinned data and 
fits a and (3 simultaneously, not separately to the size and 
period distributions. However the comparison of uncor- 
rected counts to a projected PLDF does mimic the fitting 
procedure, which applies the efficiencies to the "theory" 
(i.e. powerlaw) not the detections. 

We now explain how the powerlaw slopes are altered 
by the projection into observational space. Though the 
best fit period distribution in the right panel of Fig. 3 has 
(3 = 1.3, the curve in observational space has /3 bs — 0.3, 
smaller by about one. A decrease of precisely 2/3 is due 
to the transit efficiency. The remainder (here ~ —1/3) 
is due to the period dependence of ?7disc and depends 
on the amount of small planets in the sample. From 
Equation (4), this correction could be as large as —1.07, 
the period exponent of 7/disc at low efficiencies. Though 
not easily visible, there is modest curvature to the pro- 
jected period powerlaw due to lower discovery efficiencies 
at longer periods. 

The discovery efficiency has a much more dramatic ef- 
fect on the projected size distribution, whose curve is 
evident in the left panel of Fig. 3. There are no free 
parameters to adjust either the position or angle of the 
break, because r^isc is fixed. Above ~ 3i? ffi , the Ke- 
pler discovery efficiency is unity for all periods consid- 
ered. Thus the slope of the projected size distribution 
at large radii matches the bestfit a. The break starts at 
R < 3i?©, because rydisc drops below 90% for these sizes, 
as shown in Fig. 2. At small R, the powerlaw slope of the 
projected radial PLDF is a + 6.34, as fixed by the radial 
dependence of the discovery efficiency in Equation (4) . 

The thickness of the fits shown in Fig. 3 (and other 
figures) indicate la deviations from maximum likelihood. 
Fig. 6 shows the corresponding error ellipse in the a and 
(3 parameters (the filled black oval is for the full planet 
sample). The la errors on a and (3 reported in Table 1 
arc one-dimensional errors holding the other parameter 
fixed. As explained in Section 3, the normalization C is 
not independently varied during the fit procedure, but 
follows from the values of a, (3 and the planet- to-star 
ratio. The errors on C indicate the range of values taken 
inside the la error ellipse of a and (3. 

Our best fit size distribution has a = —1.73 ± 0.07. 
This value indicates that smaller planets are much more 
abundant, but that larger planets contain most of the 
mass. At constant planet density, a < — 3 would give 
small planets most of the mass. However larger planets 
tend to have lower density, p. Take p oc R~ - 7 as an 
example. 6 Then a < —2.3 would give small planets more 
of the mass (over the range where the density law holds). 
The conclusion that larger planets dominate the mass 
distribution holds even if the giants are quite inflated. 

While planet counts drop towards small sizes, so does 
the projection of the bestfit powerlaw, as described 
above. The good qualitative agreement between the 
counts and the curve (in the left panel of Fig. 3) means 

6 Appropriate if 1 i?0 planets have p « 5 g/cm 3 and 10 i?0 
planets have p ss 1 g/cm 3 . 
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Planet Radius [i? Earth ] Period [days] 

Fig. 3. — The histograms give uncorrected Kepler planet candidate counts around Solar type stars. Counts are binned by planet radius, R 
(left) and orbital period, P (right) and reported per star, normalized to bin width (i.e. divided by A In Ft or A In P, respectively) . Error bars 
measure the square root of planet counts per bin. The main filled histogram is for detections with S/N > 10, while the dotted histogram 
only includes detections with S/N > 20 (as reported by BKepll). Black curves show the best fit powerlaw oc R a pP . This powerlaw is 
convolved with the selection effects and marginalized over P (or R) before plotting against the planet counts. The la error ellipse of the 
maximum likelihood fit covers a = —1.73 ± 0.07 and ft = 1.22 ± 0.04 as shown by the wide red curve. The drop in planet counts below 
~ 2i?g) is primarily due to selection effects, as evidenced by the corresponding drop in the powerlaw fit when projected into observational 
space. The deficit of planets at P < 3 days is real and not explained by any selection effect. Because of this short period deficit, the period 
distribution is poorly described by a single unbroken powerlaw. 
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Fig. 4. — Similar to Fig. 3 but the planets are divided into a short period (P < 7 days, in green) and a longer period (P > 7 days, in 
orange) sample. Each sample is independently fit to a powerlaw PLDF (oc R a P@). Continuity at the period boundary is not enforced 
so that one sample does not affect the other. The best fit powerlaw indices are labelled "fast" and "slow" for the short and long period 
samples, respectively. The difference in the period distributions (/3) has extremely high significance. The radius powerlaw a is steeper for 
the longer period sample, a result with 2.3cr significance. Thus the ratio of small to big planets is higher at longer periods. 



that the drop in planet counts is mostly due to selection 
effects. The actual size distribution continues to rise to- 
wards smaller planets. 

The best fit period law, (3 — 1.33 ±0.03, indicates that 
planets are more closely packed further from the star, 
for logarithmic intervals in P or semimajor axis a. For 
/3 > 3/2 the planet density per linear interval in a would 
increase. 



However Fig. 3 shows that a single powerlaw is a quali- 
tatively poor fit to the data due to a sharp drop in planet 
counts below P ~ 3 days. Planet counts flatten at longer 
periods. As explained above, this means that the actual 
period distribution continues to rise towards longer peri- 
ods. This break in the period distribution motivates our 
division of the planet sample below. 
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Fig. 5. — The planet sample is divided into four quadrants with separations in period at 7 days and in radius at 3i?0. Each quadrant is 
separately fit to a powerlaw PLDF (oc R a pP). The result is plotted in observational space, and the values of the best fit exponents are 
shown. The differences between the subsample distributions and their significance is discussed in the text and shown in Figs. 6 and 7. 
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Fig. 6. — Error ellipses (1,2 & 3<j) for the powerlaw fits to the 
PLDFs (oc R a PP) of different planet samples. The black contours 
show the fit to the full planet sample. Fits to the four quadrants 
of the planet sample shown in Fig. 5 are also given. The differ- 
ence in the period distribution of short and long period planets is 
highly significant. The behavior of the size distributions is more 
complex and has ~ 2 — 3<r significance. For shorter periods, the 
size distribution steepens (more negative a) between the large and 
small planet samples. Longer periods show the opposite trend — 
the size distribution flattens going from bigger to smaller planets. 

4.2. Short Vs. Long Periods 

Since the period distribution deviates strongly from a 
simple powelaw, we divide the planets into short and long 
period samples. We use P — 7 days as the dividing line. 
The choice gives us comparable numbers of planets in 
each sample and also gives an adequate range of periods 
over which to measure a powerlaw slope. Planet counts 
in different cuts are summarized in the N p \ column of 
Table 1. 



Fig. 7. — Size distributions, marginalized over period. Fits are 
to a powerlaw PLDF, for both our full planet sample and the four 
quadrants shown in Figs. 5 and 6. Planets around ~ 3i?0 appear 
preferentially depleted at orbital periods below ~ 7 days. This 
deficit could be due to the loss of volatiles from lower mass giant 
planets that approach their host stars too closely. 

Fig. 4 shows the planet counts and the results of in- 
dependent powerlaw fits to the short and long period 
data. These fits are not broken powerlaws, as the val- 
ues are not required to match at the period boundary. 
Two powerlaws vastly improve the qualitative fit to the 
data. More complicated functions could more precisely 
describe the location and nature of the period break, but 
are not considered here. 

The radius distribution in the left panel includes over- 
lapping histograms for the short and long period sam- 
ples. The longer period ("slow") sample has a steeper 
radius distribution than the shorter period ("fast"). The 
difference of a s iow — «fast ~ 0.5 ± 0.2 has about 2.5a sig- 
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nificance. By eye, the difference in the distribution is 
clearly driven by the flatness of the fast sample between 
~ 3 — 15i?0. To clarify this behavior we now examine the 
differences between small and large radius planet sam- 
ples. 

4.3. Quadrants 

We subdivide the planets into four samples by splitting 
both the shorter and longer period populations (still de- 
fined relative to 7 days) into "small" and "big" samples 
relative to The results of the four independent 

powerlaw fits are plotted against planet counts in Fig. 5. 
Fig. 6 plots error ellipses to compare the slope changes 
and their statistical significance. The error ellipses of the 
subsamples are much larger than the full sample. Smaller 
number counts and the decreased leverage of the reduced 
range in R and P (over which to define a slope) both play 
a role. 

The period distributions are only modestly affected by 
this split. At both short and long periods, the best fit j3 
of the small and big samples agree within the statistical 
uncertainties. Nevertheless some qualitatively interest- 
ing features appear in the right panel of Fig. 5. The 
small planet counts begin their decline below about 7 
days (as noted by HKepll). The large planet sample re- 
mains flat towards shorter periods — with evidence of a 
peak at P w 3 days. This peak is responsible for the flat 
appearance of the overall period distribution down to 3 
days, as seen in the right panels of Figs. 3 and 4. 

The behavior of the size distributions is even more in- 
triguing. For big planets, the best fit radius law changes 
by a b ig, s iow - "big.fast = -1.31 ± 0.47, with almost 3cr 
significance. Since the discovery efficiency of Kepler is 
near unity for these large planets, no obvious systematic 
uncertainties should be at work. For small planets, the 
change in the size distributions, a sm aii, s iow - awii.fast = 
0.70 ± 0.47. While of more modest significance, this 
change is of the opposite sign. 

Fig. 7 plots the best fit size distributions for each quad- 
rant, with the full planet sample shown for comparison. 
At longer periods the size distribution flattens towards 
smaller radii, i.e. a is less negative for the small planets. 
For the shortest period planets the behavior is precisely 
opposite. The size distribution is quite shallow at large 
sizes and steepens for the small radius sample. 

The implications of this behavior are intriguing. The 
processes that form and/or migrate planets inside ~ 7 
days clearly disfavors planets of a few to several R®. A 
plausible scenario (but lacking in specifics) follows. At 
short periods, only the most massive giants — true Jo- 
vians with R > 11 R® can efficiently retain their atmo- 
spheres. Lesser giants are stripped of their atmospheres, 
including ices that sublimate from the surface. The de- 
motion of ice and gas giants to small rocky cores can si- 
multaneously explain the flatter size distribution at large 
sizes and the steeper distribution at small sizes. 

If this hypothesis is born out by further observations 
and theoretical study, it would provide yet another pil- 
lar of support for the core-accretion hypothesis and for 
migration as the source of short period planets. 

4.4. Systematic Uncertainty 

To test the quality of the fits, we consider an increased 
S/N threshold. Increasing this threshold preferentially 



removes smaller planets from the sample. In principle the 
appropriate S/N threshold is set by the analysis of the 
discovery efficiency, ?7disc- As discussed in §2.3, increasing 
the S/N threshold is a conservative way to check how the 
assumed rfaisc impacts the results. 

Fig. 3 includes dotted histograms that correspond to 
the higher S/N > 20 sample. Planet counts are notice- 
ably lower below 3i?© and across all periods. 

Table 1 includes fits to the higher S/N planet sam- 
ples. All of the qualitative trends discussed in this section 
persist in the higher S/N sample and with only mostly 
reduced significance. Indeed, the interesting trends we 
identified for small planets was a steepening of the size 
distribution at small periods and a flattening at longer 
periods. It would be rather difficult for systematic un- 
certainties to conspire to produce both trends. 

Despite this robustness, ?7disc is crucial for the statis- 
tical analysis of small planets. The final rows of Table 
1 shows that setting 77di sc = 1 (but still including the 
transit probability) gives very discrepant results. These 
discrepancies are expected due to the low discovery effi- 
ciency below 2i?®, as shown in Fig. 2. 

5. PLANET OCCURRENCE 

To derive the number of planets per star, NPPS, we 
consider our non-parametric analysis. For R > 2.0i?© 
and P < 50 days, we find /2,5o = 0.19 planets per star, 
consistent with HKepll. 

When we consider smaller planets the NPPS is roughly 
one. For R > 0.5i?©, S/N thresholds of 10, 15 and 20 
give /o.5,5o = 1-36, 1.05 and 0.72 planets per star. Since 
the S/N = 20 threshold is conservative, it seems likely 
/o.5,50 > 1- 

The likely existence of more than one planet per So- 
lar type star, especially within 50 days, is remarkable. 
However it does not imply that the Sun is exceptional 
for lacking a planet so close. For R > 1.5i?© there 
are ^2 — 3 planets per planetary system (Lissauer et al. 
2011b). Multiplicity rates will almost certainly increase 
when planets down to 0.5i?© are considered. 

It is inherently speculative to extrapolate into unob- 
served regions of parameter space. Nevertheless, all the 
powerlaw fits of the previous section have distributions 
that rise with period. To emphasize the implications of 
this rise, we calculate the extrapolated density of Earth- 
sized planets at 1 AU. The final column of Table 1 ex- 
presses the normalization to the fits as 

^" ain^lnP ^'^' (38) 

the NPPS in a logarithmic interval centered on an Earth 
radius and a year. 

The most relevant fit is to the small planet, long period 
sample, which gives C©. yr = 2.75 ±0.33. Integrating this 
powerlaw gives on average /1,36s ~ 3 Earth-like planets 
with periods below a year. If the distribution function 
of planets does not turn over shortly beyond periods of 
50 days, the implications for habitable Earths around 
Sun-like stars is truly staggering. 

6. CONCLUDING DISCUSSION 

We develop a general technique for the analysis of exo- 
planet survey data. The merits of the approach, a gener- 
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alization of the TT02 maximum likelihood analysis, in- 
clude the following: 

• Data are analyzed without binning in order to pre- 
serve the statistical significance of each detection. 

• Describing the number of planets per star (NPPS) 
allows multi-planet systems to be included. 

• Selection effects are used to calculate the number 
of expected detections. They are not used to en- 
hance (or diminish) the number of actual detec- 
tions, which affects error estimates. 

• The overall normalization is analytically removed 
from the likelihood function. With one less param- 
eter in the numerical optimization, the fit is simpler 
and possibly more robust. 

• Different surveys can be jointly analyzed, as shown 
by TT02. Even different survey methods that 
probe different regions of parameter space can be 
combined if selection effects are well characterized. 

We apply the technique to Kepler planet candidates 
released by BKepll. We focus on planets orbiting Solar 
type stars for which HKcpll has quantified the Kepler 
discovery efficiency. Fits to these efficiences (see Fig. 2) 
reveal a remarkably smooth variation with planet radius 
and orbital period that agrees with analytic predictions 
for transit surveys. This regular behavior gives us confi- 
dence to analyze detections down to 0.5i?©. 

From this analysis, the best estimate of the number 
of planets per star bigger than 0.5i? ffl and with periods 
below 50 days is / ~ 0.7 — 1.4. Uncertainty in the dis- 
covery efficiencies dominate these estimates. Since planet 
occurrence rises towards both small sizes and longer peri- 
ods, the only question is when (not if) we can claim with 
certainty that there is more than one planet per star. 

The shape of the radius and period distributions is even 
more informative than absolute number counts. The 
most notable trend is a difference between the size dis- 
tributions of shorter and longer period planets (below 
and above 7 days), plotted in Fig. 7. While the shorter 
period planets are less abundant overall, their relative 
deficit is most pronounced around <~ 3i? ffl . We interpret 



this finding as arising from the preferential evaparation 
of ice and lower mass gas giants that migrated too close 
to their host star. In this picture, the steepness of the 
size distribution of small planets at short periods arises 
from the remnant cores of these stripped giants. 

If this interpretation holds, it would add support to the 
core accretion hypothesis (Pollack et al. 1996), though 
this theory is not seriously challenged at masses below 
Mj up and possibly much higher (Kratter et al. 2010). In 
particular, cores would have to be present by the end of 
migration. This constraint applies to models where the 
cores sediment over time (Helled & Schubert 2008). 

We do not here attempt to constrain the mode of 
migration with these size trends. Possibilities include 
smooth inward migration through the disk, a theory that 
was highly developed even before the discovery of extra- 
solar planets (Lin & Papaloizou 1986; Ward 1986; Arty- 
mowicz 1993). Planet- planet scattering (Chatterjee et al. 
2008) and Kozai oscillations (Wu & Murray 2003) have 
also been proposed as mechaisms to deliver planets close 
to their hosts. Recently Wu & Lithwick (2010), proposed 
a new mechanism — secular chaos — and compared all 
the proposed mechanisms to observations of hot Jupitcrs. 
These theories and others should now be compared to 
Kepler data for low mass planets as well. 
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APPENDIX 

A. ANALYTIC SOLUTIONS FOR POWERLAW FITS 

This appendix derives an analytic solution to the best fit powerlaws of the planetary radius and orbital period. The 
transit probability is included, but all other detection efficiencies are set to unity, as appropriate for a perfect transit 
survey. Since numerical fits are more flexible, this derivation is intended to provide insight and to provide a check on 
numerical solutions (when only the transit probability is included as a selection effect). 

Consider the detection of N p \ planets in a transit survey. We fit the PLDF to the powerlaw of Equation (31) 
(equivalent to Equation (36) used in our main analysis of Kepler data). The logarithms of radius and period are 
defined as the parameters x and y as in Equation (32). We define the reference planet radius, R , and orbital period 
P so the (here rectangular) domain of the survey is centered on x = y = and all detections (indexed by i) fall within 

-x m <Xi<x m ; ~y m <yi<ym- (Al) 

The transit efficiency of Equation (1) can be written as 



1/3 

Vtv = ?7tr,oCXp ( ) ; Tftr.o = ( Q p "p2 ) ' ( A2 ) 



2y\ _ ( 3?r 



and we ignore other selection effects. 
The number of expected transits for such a perfect survey is N exp = N*CF peT f with 



sinh(ax m ) sinhjQg - 2/3)y m ] 
Fpei{ = 4 VtI , j-^ . (A3) 
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given by Equation (29) (which averages the general Equation (17) over stellar parameters). 
The best-fitting a and j3 follow from Equation (23) as 

coth(ax m ) - J- = ; coth[(/3 - 2/% m ] - ^- = . ( A4 ) 

Mm a; m (P ~ 2/3)j/ m j/ m 

where the obsevational mean (:r)obs = ^ Ej^ii an< ^ similarly for (y) bs- The relevant transcendental function, 
T(7) = coth(7) — I/7 is monotonic in 7 with T — > ±1 as 7 — >• ±00 and T(7) « 7/3 for I7 <C 1. The closed form 
solutions of Equation (A4) have several features that relate to more complete analyses: 

• The powerlaws increase monotonically with the average value of the (log of the) observables, independent of how 
the values are distributed around that mean. Higher order shape functions can describe more detailed features 
in the data. 

• The period powerlaw f3 in the PLDF is 2/3 larger than the period law describing the detections. The transit 
efficiency, r] tI causes this increase. In general any powerlaw selection effect adjusts from the observed to the 
underlying distribution this way. 

• The overall magnitude of the detection efficiency does not affect powerlaws, or shape parameters in general. 
For instance, uncertainty in r/ t r,o due to uncertainty in stellar densities does not affect the shape function. The 
overall planet occurrence of course does depend on the magnitude of detection efficiencies, via the normalization 
given by Equation (21). 

• Solutions for the best fit powerlaws decouple. In general, the solution for two shape parameters decouple if the 
shape parameters (and the physical parameters they describe) are separable in the shape function g and the 
relevant physical parameters are separable in the detection efficiencies. This separability condition is not met for 
out full treatment of the Kepler data. While the efficiency ratio rdi sc of Equation (4) is separable, the efficiency 
itself, ?7di SC = rdisc/ (1 + rdi sc ), is not. This leads to the modest covariance between a and /3 seen in Fig. 6. 
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TABLE 1 

Joint Powerlaw Fits to Kepler Planet Candidates in "Solar'' Subsample 







Planet Sample 








Fit Properties 




Tag 


Radii 


Periods 


Wpi b 


JVtot 


plancts c 


a 


P 






[«el 




[days] 




(corrected) 


per star 


(radius law) 


(period law) 




Full 


0.5 — 


20 


0.5 — 50 


562 


1.4 X 10 5 


2.38 


-1.73 ± 0.07 


1.33 ± 0.03 


23.40 ± 0.50 


Full/K2 


2 — 


20 


0.5 — 50 


372 


1.2 X 10 4 


0.21 


-1.88 ± 0.11 


1.22 ± 0.04 


20.30 ± 1.35 


(Short vs. Long Period: 


"Fast" vs. ' 


Slow") 












Fast 


0.5 — 


20 


0.5 — 7 


211 


1.0 x 10 4 


0.18 


-1.44 ± 0.11 


2.23 ± 0.12 


(15 ± 1)E2 


Slow 


0.5 — 


20 


7 — 50 


351 


1.3 x 10 5 


2.26 


-1.93 ± 0.10 


0.63 ± 0.10 


3.57 ± 0.06 


Fast/i?2 


2 — 


20 


0.5 — 7 


111 


1.3 x 10 3 


0.02 


-1.09 ± 0.17 


2.27 ± 0.18 


(9.8 ± 1.3)E2 


Slow/K2 


2 — 


20 


7 — 50 


261 


9.7 x 10 3 


0.17 


-2.31 ± 0.15 


0.54 ± 0.11 


4.71 ± 0.47 


Small 


0.5 - 


- 3 


0.5 — 50 


389 


1.3 x 10° 


2.30 


-1.52 ± 0.16 


1.43 ± 0.04 


32.05 ± 2.16 


Big 


3 — 


20 


0.5 — 50 


173 


4.5 x 10 3 


0.08 


-1.42 ± 0.16 


1.08 ± 0.06 


5.26 ± 0.32 


(Quadrants) 




















Small/Fast 


0.5 - 


- 3 


0.5 — 7 


148 


1.5 x 10 4 


0.25 


-1.93 ± 0.24 


2.35 ± 0.14 


(33 ± 3)E2 


Big/Fast 


3 — 


20 


0.5 — 7 


63 


6.8 x 10 2 


0.01 


-0.66 ± 0.24 


2.08 ± 0.22 


(1.7 ± 0.3)E2 


Small/Slow 


0.5 - 


- 3 


7 — 50 


241 


7.6 x 10 4 


1.31 


-1.23 ± 0.23 


0.68 ± 0.12 


2.75 ± 0.33 


Big/Slow 


3 — 


20 


7 — 50 


110 


3.4 x 10 3 


0.06 


-1.97 ± 0.23 


0.47 ± 0.17 


2.09 ± 0.23 



(S/N>20) 



Full/sn20 


0.5 — 


20 


0.5 — 


50 


453 


8.0 x 10 4 


1.37 


-1.44 ± 0.07 


1.30 ± 0.04 


12.54 ± 0.31 


SmFa/sn20 


0.5 - 


- 3 


0.5 - 


- 7 


105 


6.5 x 10 3 


0.11 


-1.34 ± 0.30 


2.16 ± 0.16 


(7.3 ± 0.8)E2 


BiFa/sn20 


3 — 


20 


0.5 - 


- 7 


62 


6.7 x 10 2 


0.01 


-0.51 ± 0.25 


2.07 ± 0.22 


(1.3 ± 0.2)E2 


SmSl/sn20 


0.5 - 


- 3 


7 — 


50 


178 


3.6 x 10 4 


0.62 


-0.63 ± 0.30 


0.77 ± 0.14 


1.73 ± 0.28 


BiSl/sn20 


3 — 


20 


7 — 


50 


108 


3.4 x 10 3 


0.06 


-1.88 ± 0.24 


0.50 ± 0.17 


1.98 ± 0.22 



(TRANSIT PROBABLILITY ONLY, Perfect Discovery Efficiency, r] diBC = 1) 



Full/tp 


0.5 — 


20 


0.5 — 


50 


562 


1.4 x 10 4 


0.24 


-0.18 ± 0.04 


1.04 ± 0.03 


0.64 ± 0.01 


SmFa/tp 


0.5 - 


- 3 


0.5 - 


- 7 


148 


1.5 x 10 3 


0.03 


1.19 ± 0.18 


1.94 ± 0.14 


41.49 ± 3.59 


BiFa/tp 


3 — 


20 


0.5 - 


- 7 


63 


6.8 x 10 2 


0.01 


-0.52 ± 0.25 


2.08 ± 0.22 


(1.4 ± 0.2)E2 


SmSl/tp 


0.5 - 


- 3 


7 — 


50 


241 


7.0 x 10 3 


0.12 


2.45 ± 0.18 


0.26 ± 0.12 


0.02 ± 0.00 


BiSl/tp 


3 — 


20 


7 — 


50 


110 


3.4 x 10 3 


0.06 


-1.89 ± 0.24 


0.48 ± 0.17 


1.89 ± 0.21 



a Planet radii are computed form the products of R/R* and R# given in Tables 1 and 2 of BKcpll. 

b 7V p i :Numbcr of Kepler planet candidates with S/N > 10 (or 20 where indicated) as reported in BKcpll . 

c See §5 for more accurate non-parametric fits to the planet occurrence. 



